 clear;
load('datastaticB')

phiI0=phiI;
phiX0=phiX;

for j=1:length(ss)-1
s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);
ix=ixx(j);
x=xx(j);
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j);   
end
errorode(j+1)=0;


for j=flagstar-1:-1:1
%bs(j)=bs(j+1)-dbs(j+1)*ds;
s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);
ix=ixx(j);
x=xx(j);

for kkj=1:600
db=(bs(j+1)-bs(j))/ds;

for kj=1:10
    


ixnew=(1-1/(bs(j)^(1-1/psi)*(A-ix-x)^(1/psi)/rho))/etaI;
ix=99/100*ix+ixnew/100;


end


errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j);
while abs(errorode(j))>10^(-11)
derrorode=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi-1)*(1-1/psi)*(A-ix-x)/bs(j)^2*(-1))+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*((-1)/ds/bs(j)+(bs(j+1)-bs(j))/ds/bs(j)^2*(-1));
bs(j)=bs(j)-errorode(j)/derrorode;
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j);   
end

%for kj=1:20
 xnew=(1-bs(j)/db*(1-etaI*ix))/etaX*s;
 x=9999/10000*x+xnew/10000;




end









ixx(j)=ix;
xx(j)=x;


end







for j=1:flagstar-1
s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);
ix=ixx(j);
x=xx(j);    
    
db=(bs(j+1)-bs(j))/ds;
errori(j)=(1-1/(bs(j)^(1-1/psi)*(A-ixx(j)-xx(j))^(1/psi)/rho))/etaI-ixx(j);
%errorx(j)=(1-(bs(j)-s*db)/db*(1-etaI*ixx(j)))/etaX*ss(j)-xx(j);
errorx(j)=(1-bs(j)/db*(1-etaI*ixx(j)))/etaX*ss(j)-xx(j);
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j);    
end
errori(j+1)=0;
errorx(j+1)=0;






continuedynamicB;


% subplot(3,2,1)
% plot(ss,bs,'LineWidth',1.5)
% hold on;
% 
% % plot(ss,bs0,'k:','LineWidth',1.5)
% % hold on;
% 
% plot(ss,bx,'r--','LineWidth',1.5)
% hold on;
% 
% 
% 
% 
% subplot(3,2,2)
% plot(ss,ixx,'LineWidth',1.5)
% hold on;
% 
% plot(ss,ixx0,'k:','LineWidth',1.5)
% hold on;
% 
% subplot(3,2,3)
% plot(ss,xx,'LineWidth',1.5)
% hold on;
% 
% plot(ss,xx0,'k:','LineWidth',1.5)
% hold on;
% 
% 
% subplot(3,2,4)
% plot(ss(1:flagstar),errorode(1:flagstar),'LineWidth',1.5)
% hold on;
% plot(ss(1:flagstar),errori(1:flagstar),'r--','LineWidth',1.5)
% hold on;
% plot(ss(1:flagstar),errorx(1:flagstar),'k:','LineWidth',1.5)
% hold on;
% 
% subplot(3,2,5)
% plot(ss,ixx-ixx0,'LineWidth',1.5)
% hold on;
% 
% subplot(3,2,6)
% plot(ss,bs-bs0,'LineWidth',1.5)
% hold on;

% subplot(3,2,5)
% plot(ss,dbs-dbsan,'LineWidth',1.5)
% hold on;
